Fit LFI models

Author

Max Lindmark

Published

2025-03-13

Load libraries

# Load libraries
library(tidyverse)
library(tidylog)
library(sdmTMB)
library(patchwork)
library(viridis)
library(RColorBrewer)
library(modelr)
library(ggstats)
library(ggtext)
library(ggspatial)
library(ggsidekick)
theme_set(theme_sleek())

home <- here::here()

source(paste0(home, "/R/map-plot.R"))
Reading layer `StatRec_map_Areas_Full_20170124' from data source 
  `/Users/maxlindmark/Dropbox/Max work/R/pelagic-lfi/data/shapefiles/ICES-StatRec-mapto-ICES-Areas/StatRec_map_Areas_Full_20170124.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 11074 features and 17 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -44 ymin: 36 xmax: 69 ymax: 85
Geodetic CRS:  WGS 84

TODO: knitting error

Read data and prediction grid, scale variables

# Read data

d <- read_csv(paste0(home, "/data/clean/lfi_data.csv")) |> 
  mutate(year_f = as.factor(year),
         log_depth = log(depth),
         log_depth_sc = (log_depth - mean(log_depth)) / sd(log_depth),
         log_depth_sq = log_depth_sc^2,
         year = as.integer(year)) |> 
  filter(!year == 1991)

names(d)
 [1] "year"              "haulno"            "date"             
 [4] "haul_id"           "sub_div"           "ices_rect"        
 [7] "lat"               "long"              "times"            
[10] "period"            "bdepth"            "tdepth"           
[13] "gear"              "gear_type"         "utc"              
[16] "tot_kg_hour"       "large_tot_kg_hour" "lfi"              
[19] "lon"               "X"                 "Y"                
[22] "depth1"            "depth2"            "depth"            
[25] "year_f"            "log_depth"         "log_depth_sc"     
[28] "log_depth_sq"     
t <- d |> drop_na(tdepth)
t <- d |> drop_na(utc)

d$utc
18:28:00
07:30:00
02:25:00
03:05:00
01:30:00
17:35:00
20:35:00
18:05:00
04:30:00
22:15:00
03:45:00
01:25:00
12:25:00
14:30:00
12:00:00
00:52:00
07:48:00
15:30:00
08:15:00
22:45:00
18:10:00
08:57:00
09:50:00
13:58:00
04:29:00
12:50:00
16:20:00
07:30:00
07:35:00
11:05:00
06:17:00
17:15:00
12:35:00
19:20:00
03:10:00
18:15:00
02:06:00
08:05:00
14:20:00
06:30:00
02:45:00
07:30:00
04:29:00
18:45:00
10:00:00
22:40:00
04:29:00
15:05:00
07:45:00
10:44:00
05:03:00
04:37:00
08:20:00
12:00:00
13:30:00
10:07:00
06:10:00
08:46:00
11:55:00
07:25:00
16:40:00
20:25:00
02:20:00
18:20:00
15:13:00
14:15:00
19:54:00
22:45:00
14:30:00
05:54:00
02:15:00
13:30:00
04:37:00
18:42:00
06:40:00
10:46:00
17:28:00
07:27:00
08:29:00
18:31:00
13:15:00
06:58:00
19:39:00
06:25:00
19:50:00
07:00:00
08:12:00
17:13:00
04:42:00
17:18:00
11:34:00
19:29:00
18:17:00
18:31:00
04:22:00
06:36:00
12:59:00
12:59:00
07:12:00
17:52:00
19:09:00
19:50:00
14:15:00
17:09:00
07:55:00
20:23:00
03:51:00
16:25:00
18:38:00
05:30:00
19:25:00
17:15:00
03:50:00
08:05:00
12:05:00
17:27:00
21:45:00
15:00:00
14:00:00
17:10:00
12:10:00
07:50:00
03:30:00
02:00:00
13:30:00
06:00:00
21:37:00
10:35:00
01:32:00
08:49:00
19:43:00
07:20:00
18:59:00
21:45:00
10:01:00
23:36:00
18:15:00
02:50:00
02:16:00
10:35:00
12:50:00
23:35:00
16:08:00
20:43:00
17:15:00
18:34:00
19:55:00
13:12:00
07:30:00
06:15:00
10:30:00
07:16:00
17:27:00
08:05:00
14:15:00
14:25:00
20:28:00
15:15:00
18:36:00
10:30:00
07:13:00
16:50:00
10:58:00
18:30:00
14:10:00
10:28:00
10:50:00
12:31:00
10:30:00
12:00:00
13:22:00
10:49:00
14:05:00
07:42:00
10:04:00
08:08:00
21:15:00
10:58:00
18:50:00
04:40:00
15:16:00
12:54:00
07:30:00
10:25:00
09:35:00
14:02:00
11:46:00
15:21:00
12:45:00
19:15:00
09:41:00
13:49:00
17:35:00
09:58:00
02:45:00
15:33:00
09:37:00
11:45:00
11:21:00
08:02:00
18:45:00
13:08:00
16:19:00
13:48:00
08:17:00
17:34:00
10:20:00
09:44:00
08:07:00
07:00:00
17:23:00
11:42:00
10:09:00
08:44:00
09:14:00
18:02:00
13:53:00
15:20:00
14:17:00
05:38:00
10:25:00
11:40:00
13:36:00
15:42:00
00:50:00
19:20:00
14:30:00
17:00:00
17:40:00
15:20:00
18:30:00
10:35:00
21:00:00
16:17:00
09:15:00
07:35:00
15:30:00
17:10:00
14:30:00
13:40:00
07:28:00
21:10:00
01:20:00
19:55:00
12:45:00
16:15:00
16:35:00
10:45:00
21:58:00
14:50:00
15:24:00
18:06:00
07:45:00
20:22:00
04:28:00
16:08:00
13:59:00
12:24:00
05:58:00
13:23:00
19:50:00
12:26:00
18:11:00
18:10:00
00:25:00
15:34:00
04:34:00
14:37:00
15:30:00
05:01:00
07:28:00
08:53:00
13:50:00
15:37:00
20:25:00
16:55:00
14:45:00
09:54:00
17:01:00
12:23:00
20:35:00
05:16:00
21:30:00
11:40:00
04:46:00
19:08:00
17:25:00
05:30:00
18:49:00
04:43:00
05:47:00
04:24:00
09:49:00
14:40:00
12:55:00
18:35:00
00:05:00
20:08:00
19:19:00
11:15:00
21:50:00
21:22:00
05:18:00
18:31:00
00:40:00
11:50:00
16:15:00
16:50:00
19:44:00
13:20:00
15:55:00
20:37:00
17:59:00
06:30:00
14:30:00
19:30:00
07:15:00
07:20:00
19:45:00
19:02:00
06:20:00
08:35:00
21:31:00
16:17:00
21:09:00
07:20:00
02:44:00
19:07:00
14:18:00
09:10:00
04:31:00
05:34:00
12:05:00
17:00:00
05:15:00
04:48:00
04:30:00
18:29:00
07:15:00
10:15:00
05:27:00
05:50:00
07:06:00
22:22:00
20:33:00
02:00:00
00:15:00
05:40:00
08:24:00
08:53:00
04:20:00
07:25:00
17:20:00
20:05:00
03:00:00
13:56:00
21:15:00
16:02:00
07:30:00
12:00:00
21:30:00
02:15:00
09:33:00
21:30:00
09:00:00
14:45:00
04:37:00
22:15:00
13:36:00
07:44:00
10:58:00
07:58:00
20:20:00
01:50:00
03:27:00
13:26:00
07:54:00
03:34:00
14:25:00
16:00:00
17:40:00
07:30:00
01:35:00
15:35:00
02:10:00
21:39:00
11:55:00
09:25:00
23:36:00
01:10:00
17:02:00
17:33:00
12:10:00
03:20:00
13:10:00
09:09:00
00:25:00
14:11:00
13:15:00
05:41:00
09:14:00
13:38:00
15:21:00
11:47:00
05:00:00
11:57:00
07:45:00
03:33:00
11:45:00
08:25:00
18:51:00
07:16:00
13:57:00
05:15:00
15:00:00
12:55:00
09:45:00
10:45:00
04:55:00
21:28:00
13:52:00
11:10:00
09:55:00
10:57:00
06:38:00
04:24:00
08:44:00
09:03:00
06:40:00
11:55:00
06:27:00
10:15:00
10:30:00
10:00:00
10:51:00
09:35:00
06:39:00
15:17:00
19:17:00
18:20:00
17:35:00
11:50:00
13:50:00
13:29:00
13:36:00
11:25:00
09:40:00
11:19:00
15:16:00
17:14:00
09:07:00
16:45:00
21:05:00
11:05:00
19:26:00
06:21:00
04:27:00
06:05:00
00:00:00
01:39:00
02:30:00
15:19:00
17:55:00
09:05:00
20:40:00
12:42:00
04:44:00
16:15:00
06:20:00
15:42:00
18:13:00
16:15:00
16:18:00
12:00:00
10:49:00
07:20:00
22:45:00
16:38:00
08:05:00
14:08:00
20:11:00
04:20:00
14:40:00
15:25:00
11:57:00
16:39:00
10:30:00
06:35:00
21:13:00
15:59:00
13:49:00
11:55:00
22:40:00
10:55:00
01:39:00
12:59:00
18:25:00
18:00:00
19:30:00
16:45:00
09:20:00
15:51:00
14:59:00
18:13:00
14:23:00
20:29:00
21:15:00
07:31:00
05:35:00
01:50:00
20:13:00
13:30:00
04:05:00
19:50:00
07:40:00
01:06:00
02:50:00
18:57:00
01:45:00
10:08:00
21:57:00
10:50:00
08:53:00
07:45:00
02:55:00
10:40:00
13:03:00
11:19:00
06:35:00
04:19:00
15:00:00
06:27:00
07:12:00
20:02:00
12:00:00
07:57:00
07:35:00
07:45:00
18:55:00
19:41:00
04:35:00
19:22:00
18:09:00
18:00:00
17:57:00
12:58:00
21:39:00
05:09:00
10:55:00
19:30:00
07:45:00
16:16:00
06:34:00
14:25:00
11:00:00
13:13:00
07:55:00
07:50:00
01:55:00
20:07:00
09:47:00
10:56:00
14:15:00
01:28:00
14:45:00
07:00:00
20:55:00
10:58:00
02:50:00
14:40:00
11:47:00
12:00:00
12:39:00
02:20:00
05:37:00
04:37:00
15:00:00
16:02:00
17:37:00
04:34:00
13:35:00
04:35:00
19:50:00
08:06:00
20:55:00
04:48:00
15:56:00
03:48:00
09:48:00
05:16:00
10:29:00
04:10:00
09:56:00
05:12:00
13:30:00
07:35:00
18:15:00
05:30:00
07:46:00
11:22:00
17:18:00
16:03:00
02:09:00
23:25:00
18:18:00
14:38:00
09:04:00
08:42:00
19:27:00
08:33:00
14:03:00
07:50:00
02:25:00
07:37:00
06:19:00
08:05:00
20:11:00
06:59:00
16:40:00
21:20:00
20:28:00
14:27:00
20:57:00
13:33:00
19:03:00
14:20:00
07:58:00
18:10:00
06:18:00
16:10:00
19:29:00
03:15:00
23:20:00
20:00:00
21:20:00
04:29:00
20:29:00
03:18:00
12:40:00
22:19:00
18:19:00
12:04:00
04:22:00
03:25:00
17:26:00
21:09:00
09:45:00
01:35:00
06:25:00
19:54:00
12:05:00
07:25:00
04:50:00
17:13:00
16:35:00
21:15:00
13:06:00
16:22:00
08:43:00
06:18:00
06:19:00
12:05:00
01:10:00
06:33:00
00:55:00
15:54:00
11:59:00
04:35:00
16:10:00
02:05:00
19:26:00
23:40:00
11:39:00
00:15:00
14:15:00
23:00:00
18:56:00
17:52:00
08:47:00
09:10:00
17:11:00
13:13:00
16:35:00
17:05:00
04:21:00
18:30:00
17:45:00
18:29:00
13:11:00
15:14:00
21:15:00
20:55:00
13:30:00
03:28:00
22:12:00
12:38:00
14:18:00
20:25:00
17:21:00
20:26:00
13:04:00
12:00:00
21:35:00
06:48:00
18:20:00
18:59:00
20:32:00
23:25:00
16:58:00
17:50:00
14:34:00
15:41:00
14:25:00
19:00:00
07:25:00
04:56:00
10:49:00
18:35:00
07:35:00
12:35:00
17:20:00
13:35:00
05:46:00
13:10:00
10:57:00
16:25:00
19:53:00
04:43:00
14:40:00
04:38:00
08:15:00
19:20:00
05:26:00
20:40:00
00:20:00
19:57:00
09:29:00
13:36:00
08:26:00
19:57:00
04:39:00
08:30:00
21:44:00
19:58:00
20:10:00
07:41:00
17:25:00
18:54:00
08:00:00
15:23:00
08:53:00
00:40:00
11:40:00
18:54:00
04:28:00
18:30:00
04:15:00
15:40:00
04:39:00
12:25:00
07:35:00
07:12:00
15:30:00
15:00:00
05:12:00
11:55:00
03:10:00
10:47:00
22:45:00
15:57:00
16:49:00
13:09:00
15:51:00
12:45:00
17:15:00
15:48:00
17:25:00
19:55:00
10:07:00
11:55:00
15:40:00
20:09:00
07:20:00
17:57:00
06:56:00
22:55:00
18:25:00
12:35:00
13:45:00
19:19:00
03:29:00
20:26:00
01:10:00
12:52:00
16:42:00
13:10:00
13:38:00
13:00:00
14:55:00
21:21:00
11:44:00
04:27:00
12:15:00
04:36:00
07:30:00
20:16:00
07:50:00
02:45:00
04:29:00
01:45:00
15:23:00
15:45:00
15:21:00
14:35:00
11:09:00
04:26:00
10:15:00
12:00:00
02:25:00
20:13:00
17:49:00
17:16:00
09:20:00
17:04:00
20:35:00
09:59:00
17:49:00
16:27:00
13:16:00
12:05:00
11:00:00
07:55:00
12:35:00
22:50:00
17:26:00
13:00:00
12:25:00
11:25:00
04:56:00
00:00:00
18:16:00
18:57:00
12:05:00
13:22:00
06:31:00
15:22:00
03:34:00
07:00:00
11:02:00
16:35:00
17:01:00
16:30:00
20:39:00
02:15:00
19:12:00
23:27:00
18:04:00
19:45:00
04:30:00
05:09:00
22:15:00
22:50:00
19:55:00
08:40:00
07:32:00
06:20:00
17:49:00
20:39:00
18:29:00
09:44:00
06:45:00
20:17:00
09:05:00
21:35:00
15:20:00
05:34:00
16:00:00
16:37:00
21:06:00
20:55:00
20:14:00
01:40:00
04:21:00
04:34:00
07:21:00
16:24:00
07:19:00
04:37:00
19:37:00
17:30:00
02:33:00
04:31:00
21:11:00
08:30:00
18:26:00
14:25:00
19:30:00
13:06:00
15:56:00
18:49:00
08:05:00
02:50:00
14:05:00
12:50:00
12:55:00
18:30:00
05:29:00
08:40:00
22:40:00
16:45:00
19:32:00
14:40:00
09:11:00
04:28:00
11:13:00
12:42:00
13:48:00
07:50:00
15:10:00
14:01:00
15:26:00
14:35:00
16:20:00
15:11:00
12:20:00
19:05:00
09:32:00
06:30:00
08:10:00
03:40:00
07:02:00
06:14:00
11:29:00
20:10:00
21:59:00
20:46:00
11:55:00
10:55:00
16:20:00
05:09:00
08:27:00
03:48:00
13:55:00
18:33:00
07:53:00
20:50:00
02:20:00
20:15:00
09:15:00
06:20:00
07:30:00
07:45:00
09:37:00
19:27:00
17:10:00
16:50:00
13:08:00
22:45:00
07:21:00
20:15:00
15:45:00
05:45:00
21:20:00
20:15:00
14:25:00
04:37:00
07:30:00
11:00:00
20:33:00
16:22:00
03:40:00
06:22:00
16:25:00
21:05:00
13:45:00
05:04:00
19:20:00
10:11:00
04:39:00
17:05:00
11:05:00
19:24:00
02:24:00
17:35:00
17:40:00
09:24:00
21:20:00
10:40:00
01:50:00
17:50:00
04:25:00
18:36:00
12:10:00
08:05:00
21:47:00
14:30:00
08:59:00
15:22:00
08:30:00
14:14:00
07:30:00
07:30:00
20:20:00
23:41:00
20:41:00
09:50:00
19:45:00
08:05:00
04:31:00
01:20:00
00:25:00
06:31:00
12:00:00
04:28:00
02:30:00
22:25:00
12:25:00
18:15:00
22:29:00
04:46:00
11:54:00
00:02:00
15:37:00
02:00:00
21:56:00
01:50:00
07:34:00
16:42:00
20:17:00
21:35:00
03:30:00
06:59:00
08:19:00
06:52:00
13:10:00
06:25:00
15:50:00
14:18:00
10:15:00
11:33:00
09:52:00
23:27:00
04:15:00
08:36:00
07:30:00
18:13:00
04:25:00
10:40:00
04:44:00
06:03:00
23:35:00
12:00:00
13:35:00
13:36:00
19:25:00
06:30:00
10:44:00
21:36:00
10:45:00
08:42:00
07:33:00
12:25:00
04:36:00
21:21:00
12:34:00
06:49:00
05:25:00
02:45:00
16:51:00
03:20:00
05:05:00
12:59:00
12:50:00
16:45:00
20:25:00
05:00:00
14:50:00
12:53:00
08:32:00
22:25:00
04:47:00
04:45:00
09:39:00
18:21:00
16:40:00
10:24:00
02:50:00
06:03:00
06:46:00
09:17:00
17:20:00
15:12:00
04:27:00
02:10:00
13:36:00
14:16:00
16:29:00
13:20:00
04:22:00
17:30:00
00:20:00
10:53:00
05:20:00
08:43:00
06:25:00
08:06:00
10:45:00
07:03:00
18:35:00
19:22:00
11:34:00
21:11:00
08:37:00
04:25:00
20:06:00
17:21:00
17:38:00
17:00:00
05:56:00
17:19:00
11:24:00
11:09:00
04:21:00
13:45:00
10:04:00
13:43:00
17:10:00
03:29:00
14:58:00
15:20:00
04:33:00
06:31:00
06:01:00
14:02:00
14:22:00
09:54:00
14:04:00
14:00:00
04:34:00
19:29:00
09:03:00
06:55:00
11:04:00
20:20:00
20:23:00
20:38:00
15:05:00
09:22:00
08:00:00
09:39:00
20:41:00
20:20:00
12:07:00
12:16:00
06:28:00
06:30:00
16:40:00
11:20:00
16:51:00
15:46:00
12:08:00
14:19:00
11:55:00
07:35:00
07:20:00
08:15:00
20:02:00
12:31:00
23:20:00
21:25:00
13:24:00
17:48:00
12:05:00
14:20:00
07:25:00
00:02:00
18:51:00
12:11:00
17:00:00
10:05:00
21:45:00
09:38:00
06:43:00
06:47:00
17:19:00
15:29:00
14:35:00
11:48:00
08:31:00
04:36:00
13:06:00
15:11:00
04:23:00
07:00:00
18:42:00
18:47:00
06:30:00
15:42:00
12:05:00
06:37:00
03:26:00
17:28:00
14:23:00
04:25:00
11:10:00
07:14:00
11:05:00
07:30:00
08:17:00
19:19:00
21:11:00
08:48:00
14:20:00
21:47:00
04:48:00
10:20:00
04:39:00
16:26:00
14:30:00
16:15:00
19:42:00
18:01:00
20:20:00
23:00:00
12:33:00
11:15:00
05:12:00
19:15:00
06:35:00
08:22:00
20:20:00
14:05:00
22:15:00
04:31:00
07:15:00
18:28:00
21:00:00
13:55:00
05:50:00
17:55:00
01:55:00
14:14:00
11:50:00
08:00:00
00:55:00
07:20:00
21:17:00
11:55:00
07:05:00
11:07:00
15:45:00
07:50:00
15:15:00
22:15:00
12:12:00
15:55:00
18:24:00
09:25:00
15:24:00
17:15:00
17:30:00
11:45:00
03:28:00
14:05:00
02:35:00
19:50:00
06:30:00
16:36:00
13:45:00
10:14:00
08:15:00
20:50:00
23:50:00
23:15:00
20:37:00
13:30:00
21:10:00
14:30:00
17:45:00
09:15:00
04:33:00
08:55:00
06:40:00
17:31:00
16:48:00
22:03:00
15:23:00
18:58:00
16:26:00
20:48:00
17:50:00
19:09:00
13:00:00
15:24:00
17:15:00
13:02:00
12:09:00
11:18:00
15:16:00
09:02:00
21:25:00
07:07:00
11:45:00
04:45:00
05:20:00
07:45:00
01:47:00
03:51:00
20:19:00
19:52:00
03:38:00
10:51:00
18:23:00
04:48:00
09:46:00
07:58:00
16:54:00
13:28:00
09:06:00
15:12:00
06:30:00
07:05:00
12:41:00
16:27:00
09:17:00
15:57:00
13:03:00
04:44:00
11:45:00
15:35:00
08:56:00
21:56:00
19:48:00
20:02:00
13:10:00
10:49:00
06:42:00
04:34:00
16:25:00
17:31:00
13:27:00
04:28:00
04:53:00
17:39:00
09:50:00
21:12:00
07:15:00
06:35:00
06:43:00
17:47:00
11:02:00
13:55:00
07:13:00
13:58:00
19:15:00
09:45:00
13:05:00
04:44:00
05:35:00
17:42:00
20:36:00
04:59:00
09:47:00
14:29:00
21:14:00
01:38:00
17:02:00
15:10:00
06:44:00
15:13:00
10:01:00
12:30:00
14:00:00
08:47:00
08:41:00
10:25:00
19:40:00
17:05:00
12:15:00
04:49:00
17:42:00
07:11:00
19:25:00
14:04:00
08:20:00
07:25:00
20:13:00
08:59:00
13:45:00
13:52:00
15:33:00
11:07:00
11:40:00
04:34:00
12:58:00
11:55:00
02:25:00
21:56:00
18:07:00
19:24:00
13:33:00
15:32:00
06:25:00
10:54:00
11:26:00
15:39:00
16:46:00
14:26:00
09:12:00
05:07:00
12:30:00
08:57:00
17:58:00
08:46:00
11:38:00
05:05:00
07:51:00
21:10:00
04:37:00
21:28:00
15:51:00
08:35:00
06:22:00
18:24:00
04:07:00
01:12:00
20:15:00
15:18:00
14:01:00
08:56:00
19:36:00
10:50:00
15:24:00
09:24:00
18:45:00
11:15:00
13:35:00
07:15:00
06:20:00
12:04:00
02:31:00
11:45:00
07:59:00
15:27:00
11:13:00
04:33:00
06:30:00
08:58:00
13:34:00
19:15:00
15:21:00
12:59:00
04:53:00
11:50:00
04:38:00
08:50:00
17:05:00
14:53:00
07:03:00
11:21:00
05:19:00
10:34:00
07:15:00
06:30:00
07:31:00
17:24:00
04:31:00
07:56:00
10:40:00
05:03:00
15:00:00
08:44:00
19:01:00
20:50:00
08:46:00
11:03:00
11:28:00
18:31:00
06:45:00
15:57:00
04:22:00
20:29:00
18:15:00
15:25:00
13:43:00
19:10:00
13:23:00
23:37:00
16:21:00
13:48:00
04:36:00
04:39:00
12:16:00
05:41:00
07:10:00
06:39:00
10:45:00
17:05:00
05:49:00
19:00:00
16:17:00
21:39:00
00:29:00
07:02:00
12:04:00
14:32:00
16:31:00
17:43:00
23:22:00
02:32:00
16:32:00
19:56:00
23:36:00
05:06:00
14:18:00
16:58:00
22:02:00
21:48:00
01:25:00
04:44:00
08:19:00
11:53:00
02:36:00
07:04:00
11:43:00
15:34:00
23:43:00
06:07:00
12:09:00
16:19:00
20:15:00
20:44:00
00:38:00
05:03:00
08:50:00
13:06:00
18:03:00
22:22:00
02:02:00
06:17:00
12:02:00
01:27:00
16:47:00
20:32:00
01:18:00
04:38:00
15:16:00
19:03:00
22:38:00
06:30:00
15:27:00
18:44:00
07:25:00
23:40:00
05:23:00
12:37:00
16:50:00
21:13:00
00:58:00
19:01:00
00:48:00
05:08:00
10:41:00
12:12:00
16:01:00
19:48:00
02:59:00
07:04:00
16:40:00
14:51:00
20:11:00
00:49:00
06:24:00
11:21:00
18:34:00
16:47:00
23:47:00
04:05:00
08:48:00
12:56:00
18:48:00
23:15:00
02:51:00
07:47:00
12:37:00
23:31:00
18:23:00
23:34:00
03:52:00
08:54:00
14:24:00
19:35:00
00:08:00
07:22:00
10:18:00
17:24:00
03:50:00
23:08:00
04:15:00
07:07:00
12:20:00
17:28:00
20:10:00
01:06:00
05:17:00
09:19:00
13:41:00
08:03:00
17:48:00
23:45:00
04:47:00
07:56:00
15:00:00
20:45:00
01:31:00
05:04:00
11:05:00
13:48:00
18:06:00
00:13:00
03:45:00
09:29:00
14:02:00
19:57:00
18:16:00
23:28:00
03:27:00
07:56:00
12:25:00
23:02:00
17:16:00
21:38:00
03:23:00
07:53:00
12:09:00
17:58:00
23:39:00
04:15:00
10:47:00
16:58:00
04:24:00
20:53:00
10:41:00
14:23:00
18:09:00
22:58:00
04:18:00
09:43:00
15:15:00
20:11:00
01:41:00
09:24:00
07:35:00
11:17:00
16:51:00
20:01:00
00:50:00
06:35:00
11:35:00
16:24:00
19:39:00
02:18:00
14:14:00
08:39:00
15:10:00
19:17:00
21:09:00
01:15:00
05:06:00
08:31:00
11:53:00
02:47:00
17:55:00
23:09:00
02:42:00
06:45:00
11:36:00
16:24:00
20:27:00
00:38:00
03:47:00
08:48:00
07:05:00
15:32:00
19:57:00
00:40:00
07:05:00
12:29:00
17:35:00
22:47:00
02:08:00
06:15:00
10:51:00
12:00:00
16:36:00
20:53:00
01:54:00
07:09:00
11:58:00
18:17:00
23:35:00
02:35:00
06:34:00
13:27:00
16:43:00
18:03:00
22:52:00
03:03:00
07:44:00
12:04:00
19:24:00
23:52:00
03:36:00
08:33:00
13:22:00
10:22:00
14:49:00
19:03:00
00:19:00
13:48:00
14:52:00
18:05:00
02:07:00
07:21:00
11:14:00
18:16:00
23:12:00
05:00:00
10:45:00
16:32:00
23:39:00
18:18:00
03:52:00
12:33:00
14:34:00
23:21:00
04:58:00
10:49:00
16:55:00
22:47:00
03:15:00
07:54:00
22:13:00
13:30:00
19:32:00
00:34:00
06:05:00
12:24:00
17:52:00
22:49:00
03:02:00
08:00:00
10:45:00
05:03:00
16:05:00
22:51:00
02:22:00
08:35:00
12:50:00
17:40:00
00:44:00
05:59:00
10:29:00
17:26:00
23:13:00
05:15:00
09:37:00
14:53:00
19:38:00
02:14:00
06:32:00
15:34:00
12:33:00
20:35:00
01:37:00
05:06:00
10:14:00
14:11:00
17:18:00
00:38:00
04:31:00
10:18:00
16:16:00
17:02:00
00:32:00
04:58:00
10:25:00
15:19:00
23:02:00
03:14:00
07:40:00
15:08:00
19:22:00
01:40:00
00:33:00
06:31:00
10:46:00
16:47:00
22:43:00
03:37:00
10:33:00
14:28:00
18:49:00
00:57:00
04:29:00
04:56:00
10:20:00
14:33:00
19:18:00
03:49:00
08:16:00
13:00:00
18:27:00
21:53:00
23:21:00
03:40:00
08:13:00
14:05:00
18:03:00
00:51:00
04:13:00
09:27:00
14:32:00
18:22:00
03:21:00
00:55:00
05:07:00
10:34:00
05:36:00
08:16:00
14:18:00
18:43:00
01:45:00
06:22:00
10:40:00
07:26:00
16:33:00
20:44:00
01:48:00
05:57:00
12:07:00
15:48:00
21:38:00
02:49:00
06:42:00
09:15:00
13:21:00
15:31:00
18:23:00
20:57:00
02:37:00
06:24:00
11:33:00
16:41:00
19:25:00
23:44:00
07:40:00
13:19:00
00:52:00
05:18:00
11:49:00
16:44:00
22:52:00
05:42:00
12:11:00
16:55:00
22:18:00
04:03:00
08:08:00
15:28:00
20:01:00
02:34:00
06:22:00
12:57:00
18:39:00
01:07:00
05:24:00
10:37:00
14:54:00
21:21:00
01:28:00
05:21:00
21:10:00
12:14:00
16:17:00
20:08:00
23:40:00
03:16:00
08:28:00
14:28:00
19:50:00
05:14:00
12:08:00
02:31:00
16:30:00
22:22:00
05:00:00
09:24:00
16:28:00
07:02:00
10:28:00
14:14:00
14:02:00
19:19:00
23:20:00
05:11:00
13:41:00
19:51:00
00:00:00
07:51:00
14:05:00
20:54:00
02:38:00
08:28:00
15:14:00
21:53:00
02:30:00
08:01:00
15:01:00
04:58:00
08:25:00
14:28:00
20:51:00
00:55:00
06:35:00
13:34:00
19:31:00
02:17:00
08:47:00
15:39:00
20:36:00
03:30:00
09:25:00
17:13:00
21:34:00
00:55:00
07:41:00
12:51:00
20:39:00
02:06:00
08:02:00
11:29:00
20:17:00
00:27:00
06:55:00
12:48:00
16:42:00
20:17:00
03:08:00
08:13:00
13:23:00
17:24:00
00:19:00
06:48:00
13:04:00
19:20:00
00:42:00
06:41:00
10:40:00
17:02:00
23:33:00
01:50:00
06:35:00
12:19:00
18:40:00
23:27:00
06:21:00
11:01:00
17:27:00
22:08:00
02:52:00
08:05:00
15:28:00
19:13:00
00:59:00
07:06:00
14:36:00
19:18:00
02:43:00
07:35:00
12:51:00
16:05:00
23:30:00
03:34:00
07:22:00
15:22:00
20:25:00
01:05:00
05:30:00
11:49:00
16:18:00
22:41:00
02:37:00
07:12:00
13:34:00
23:41:00
04:25:00
09:26:00
14:32:00
18:24:00
22:30:00
11:13:00
16:12:00
19:56:00
00:33:00
06:17:00
10:55:00
16:26:00
23:27:00
03:57:00
08:15:00
13:17:00
05:01:00
08:00:00
13:22:00
19:45:00
00:49:00
05:21:00
09:38:00
14:23:00
19:30:00
00:27:00
03:48:00
07:51:00
15:18:00
00:58:00
05:13:00
09:31:00
15:25:00
21:49:00
04:14:00
08:53:00
14:31:00
18:36:00
00:16:00
04:48:00
10:25:00
17:26:00
21:41:00
02:49:00
08:24:00
11:06:00
15:30:00
22:54:00
01:22:00
06:23:00
13:09:00
20:17:00
04:57:00
23:37:00
06:10:00
11:56:00
19:11:00
23:33:00
04:21:00
12:42:00
16:29:00
21:27:00
00:45:00
05:06:00
11:21:00
22:15:00
03:05:00
08:46:00
13:02:00
08:00:00
18:15:00
00:43:00
07:35:00
12:19:00
16:43:00
21:12:00
03:18:00
08:59:00
14:36:00
21:15:00
01:29:00
05:45:00
12:31:00
22:38:00
01:29:00
09:21:00
16:47:00
20:55:00
00:50:00
04:59:00
10:45:00
15:53:00
04:01:00
06:50:00
10:22:00
15:35:00
05:26:00
11:41:00
16:35:00
21:33:00
04:59:00
pred_grid <- read_csv(paste0(home, "/data/clean/pred_grid.csv")) |> 
  drop_na(depth) |> 
  mutate(depth = ifelse(depth < 1, 1, depth)) |> 
  mutate(year_f = as.factor(year),
         log_depth = log(depth),
         log_depth_sc = (log_depth - mean(d$log_depth)) / sd(d$log_depth),
         log_depth_sq = log_depth_sc^2,
         year = as.integer(year)) |> 
  filter(!year == 1991) |> 
  mutate(times = mean(d$times))

Plot data

plot_map_fc + 
  geom_point(data = d,
             aes(X*1000, Y*1000, color = gear_type), size = 0.8) +
  facet_wrap(~year, ncol = 9)
Warning: `guide_colourbar()` needs continuous scales.
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).

ggsave(paste0(home, paste("/figures/map.pdf")), width = 17, height = 19, units = "cm")
Warning: `guide_colourbar()` needs continuous scales.
Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).

Fit models

mesh <- make_mesh(d, xy_cols = c("X", "Y"), cutoff = 20)

p_mesh <- ggplot() +
  inlabru::gg(mesh$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = d, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) +
  labs(x = "Easting (km)", y = "Northing (km)")
Loading required namespace: INLA
p_mesh

fit1 <- sdmTMB(lfi ~ year_f + s(times, bs = "cc") + log_depth_sc + log_depth_sq,
               data = d,
               mesh = mesh,
               family = delta_gamma(),
               spatiotemporal = "off",
               spatial = "on",
               time = "year")

fit2 <- sdmTMB(lfi ~ 1 + s(times, bs = "cc") + log_depth_sc + log_depth_sq,
               data = d,
               mesh = mesh,
               time_varying = ~ 1,
               time_varying_type = "rw0",
               family = delta_gamma(),
               extra_time = c(1991, 1993, 1997), 
               spatiotemporal = "off",
               spatial = "on",
               time = "year")

fit3 <- sdmTMB(lfi ~ 1 + s(times, bs = "cc") + log_depth_sc,
               data = d,
               mesh = mesh,
               time_varying = ~ 1,
               time_varying_type = "rw0",
               family = delta_gamma(),
               extra_time = c(1991, 1993, 1997), 
               spatiotemporal = "iid",
               spatial = "on",
               time = "year")

AIC(fit1, fit2, fit3)
     df       AIC
fit1 95 -2273.665
fit2 15 -2215.729
fit3 15 -2303.537

Check residuals

# Residuals
simulate(fit3, nsim = 500, type = "mle-mvn") |>
  dharma_residuals(fit3, plot = FALSE) |>
  ggplot(aes(observed, expected)) +
  geom_point(color = "grey30", shape = 21, size = 0.5) +
  geom_abline(col = "tomato3", linewidth = 0.6) +
  theme(aspect.ratio = 1) +
  labs(x = "Observed", y = "Expected")
Simulating ■■■                                6% | ETA: 16s
Simulating ■■■■■■■                           20% | ETA: 14s
Simulating ■■■■■■■■■■■■                      37% | ETA: 11s
Simulating ■■■■■■■■■■■■■■■■■                 54% | ETA:  8s
Simulating ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA:  5s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      88% | ETA:  2s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

ggsave(paste0(home, "/figures/residuals.pdf"), width = 11, height = 11, units = "cm")

Calculate indices

# Get index
p <- predict(fit3, newdata = pred_grid, return_tmb_object = TRUE)

ind <- get_index(
  p, 
  bias_correct = FALSE, 
  area = 1 / nrow(pred_grid |> filter(year == max(year)))
  )
filter: removed 1,051,240 rows (98%), 25,640 rows remaining
Bias correction is turned off.
It is recommended to turn this on for final inference.
# Fit trend model
p1 <- ind |> 
  mutate(period = ifelse(year > 2017, "New", "Old")) |> 
  ggplot(aes(year, est, color = period, fill = period)) + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0) +
  guides(color = "none", fill = "none") +
  scale_color_brewer(palette = "Set1") +
  geom_point(size = 2) +
  labs(y = "Predicted average LFI", x = "Year") + 
  theme(axis.title.x = element_blank()) +
  annotate("text", label = "(a)", x = Inf, y = Inf, vjust = 2, hjust = 1.5,
           color = "gray40", size = 4)
mutate: new variable 'period' (character) with 2 unique values and 0% NA
options(scipen=999)

# Trend lines
ind2 <- ind |> 
  filter(year > 2017) |> 
  mutate(cv = se / est,
         weight = 1/cv,
         weights_cs = weight / mean(weight),
         year_sc = (year - mean(year)) / sd(year),
         est_sc = est)
filter: removed 36 rows (86%), 6 rows remaining
mutate: new variable 'cv' (double) with 6 unique values and 0% NA
        new variable 'weight' (double) with 6 unique values and 0% NA
        new variable 'weights_cs' (double) with 6 unique values and 0% NA
        new variable 'year_sc' (double) with 6 unique values and 0% NA
        new variable 'est_sc' (double) with 6 unique values and 0% NA
mt <- sdmTMB(est_sc ~ year_sc,
             family = Gamma(link = "log"),
             spatial = "off", 
             #weights = weights$weights_cs,
             data = ind2)

mt
Model fit by ML ['sdmTMB']
Formula: est_sc ~ year_sc
Mesh: NULL (isotropic covariance)
Data: ind2
Family: Gamma(link = 'log')
 
Conditional model:
            coef.est coef.se
(Intercept)    -8.74    0.14
year_sc        -0.79    0.15

Dispersion parameter: 8.21
ML criterion at convergence: -50.502

See ?tidy.sdmTMB to extract these values as a data frame.
sanity(mt)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
tidy(mt)
# A tibble: 2 × 5
  term        estimate std.error conf.low conf.high
  <chr>          <dbl>     <dbl>    <dbl>     <dbl>
1 (Intercept)   -8.74      0.142    -9.02    -8.46 
2 year_sc       -0.794     0.152    -1.09    -0.495
tidy(mt, exponentiate = TRUE)
# A tibble: 2 × 4
  term        estimate conf.low conf.high
  <chr>          <dbl>    <dbl>     <dbl>
1 (Intercept) 0.000160 0.000121  0.000211
2 year_sc     0.452    0.336     0.609   
nd <- tibble(year = seq(min(ind2$year), max(ind2$year), length.out = 50)) |> 
  mutate(year_sc = (year - mean(year)) / sd(year))
mutate: new variable 'year_sc' (double) with 50 unique values and 0% NA
pp <- predict(mt, newdata = nd, se_fit = TRUE) |> 
  mutate(est = est) |> 
  mutate(lwr = exp(est - est_se*1.96),
         upr = exp(est + est_se*1.96),
         est = exp(est))
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: no changes

mutate: changed 50 values (100%) of 'est' (0 new NA)

        new variable 'lwr' (double) with 50 unique values and 0% NA

        new variable 'upr' (double) with 50 unique values and 0% NA
p2 <- ind |> 
  filter(year > 2017) |> 
  mutate(period = ifelse(year > 2017, "New", "Old")) |> 
  ggplot(aes(year, est, color = period, fill = period)) + 
  geom_line(data = pp, aes(year, est), inherit.aes = FALSE,
            color = brewer.pal(n = 3, name = "Set1")[1], 
            linetype = "longdash") +
  geom_ribbon(data = pp, aes(year, est, ymin = lwr, ymax = upr), inherit.aes = FALSE,
              alpha = 0.2) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0) +
  geom_point() +
  scale_color_brewer(palette = "Set1", direction = -1) +
  scale_fill_brewer(palette = "Set1", direction = -1) +
  guides(color = "none", fill = "none") +
  labs(y = "Predicted average LFI", x = "Year") +
  annotate("text", label = "(b)", x = Inf, y = Inf, vjust = 2, hjust = 1.5,
           color = "gray40", size = 4)
filter: removed 36 rows (86%), 6 rows remaining
mutate: new variable 'period' (character) with one unique value and 0% NA
(p1 / p2) + plot_layout(axes = "collect")

ggsave(paste0(home, "/figures/lfi_index.pdf"), width = 14, height = 17, units = "cm")

Make map plots

# Plot map predictions
plot_map +
  geom_raster(data = p$data |> filter(year == max(year)),
              aes(X*1000, Y*1000, fill = omega_s1)) +
  scale_fill_gradient2() + 
  geom_sf()
filter: removed 1,051,240 rows (98%), 25,640 rows remaining
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_raster()`).

plot_map + 
  geom_raster(data = p$data |> filter(year == max(year)),
              aes(X*1000, Y*1000, fill = omega_s2)) +
  scale_fill_gradient2() + 
  geom_sf()
filter: removed 1,051,240 rows (98%), 25,640 rows remaining
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_raster()`).

plot_map + 
  geom_raster(data = p$data, aes(X*1000, Y*1000, fill = epsilon_st1)) +
  scale_fill_gradient2() +  
  facet_wrap(~year, ncol = 8) +
  geom_sf()
Warning: Removed 42 rows containing missing values or values outside the scale range
(`geom_raster()`).

plot_map + 
  geom_raster(data = p$data, aes(X*1000, Y*1000, fill = epsilon_st2)) +
  scale_fill_gradient2() +  
  facet_wrap(~year, ncol = 8) +
  geom_sf()
Warning: Removed 42 rows containing missing values or values outside the scale range
(`geom_raster()`).

Sensitivity

Run the index only below 59.6 latitude

ggplot(d, aes(lon, lat)) + 
  geom_point() + 
  geom_hline(yintercept = 59.6)

# Get index
p2 <- predict(fit3,
              newdata = pred_grid |> filter(lat < 59.6),
              return_tmb_object = TRUE)
filter: removed 331,128 rows (31%), 745,752 rows remaining
ind2 <- get_index(
  p2, 
  bias_correct = FALSE, 
  area = 1 / nrow(pred_grid |>
                    filter(lat < 59.6 & year == max(year)))
  )
filter: removed 1,059,124 rows (98%), 17,756 rows remaining
Bias correction is turned off.
It is recommended to turn this on for final inference.
bind_rows(
  ind |> mutate(prediction = "full"),
  ind2 |> mutate(prediction = "subset")
  ) |>
  ggplot(aes(year, est, color = prediction, fill = prediction)) + 
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3, color = NA) +
  geom_line() + 
  scale_y_continuous(trans = "sqrt") +
  labs(y = "Predicted average LFI", x = "Year")
mutate: new variable 'prediction' (character) with one unique value and 0% NA
mutate: new variable 'prediction' (character) with one unique value and 0% NA

Fit with trawl depth

fit3b <- sdmTMB(lfi ~ 1 + s(times, bs = "cc") + log_depth_sc + tdepth,
                data = d,
                mesh = mesh,
                time_varying = ~ 1,
                time_varying_type = "rw0",
                family = delta_gamma(),
                extra_time = c(1991, 1993, 1997), 
                spatiotemporal = "iid",
                spatial = "on",
                time = "year")

sanity(fit3b)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large
# Get index
p3 <- predict(fit3b,
              newdata = pred_grid |> mutate(tdepth = mean(d$tdepth, na.rm = TRUE)),
              return_tmb_object = TRUE)
mutate: new variable 'tdepth' (double) with one unique value and 0% NA
ind3 <- get_index(
  p3, 
  bias_correct = FALSE, 
  area = 1 / nrow(pred_grid |>
                    filter(year == max(year)))
  )
filter: removed 1,051,240 rows (98%), 25,640 rows remaining
Bias correction is turned off.
It is recommended to turn this on for final inference.
bind_rows(
  ind |> mutate(prediction = "full"),
  ind2 |> mutate(prediction = "subset"),
  ind3 |> mutate(prediction = "tdepth")
  ) |>
  ggplot(aes(year, est, color = prediction, fill = prediction)) + 
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3, color = NA) +
  geom_line() + 
  scale_y_continuous(trans = "sqrt") +
  scale_color_brewer(palette = "Set2") +
  scale_fill_brewer(palette = "Set2") +
  labs(y = "Predicted average LFI", x = "Year")
mutate: new variable 'prediction' (character) with one unique value and 0% NA
mutate: new variable 'prediction' (character) with one unique value and 0% NA
mutate: new variable 'prediction' (character) with one unique value and 0% NA

ggsave(paste0(home, "/figures/lfi_index_sensi.pdf"), width = 14, height = 11, units = "cm")